Lévy (levy) distribution#
The Lévy distribution is a continuous distribution on a positive half-line (shifted by a location parameter) with a power-law tail.
It is the canonical one-sided \(\frac{1}{2}\)-stable law and appears naturally as:
the first-passage time of Brownian motion to a fixed level,
a heavy-tailed model for waiting times / durations with rare but enormous values,
the step-length distribution behind Lévy flights (anomalous diffusion).
A key modeling implication: the mean and variance are infinite, so moment-based intuition (sample mean, sample variance, CLT-style standard errors) is unreliable.
Learning goals#
Know the PDF and CDF (and how the CDF is written using
erfc).Understand why the mean/variance do not exist (tail behavior).
Derive and implement a NumPy-only sampler using a normal-to-Lévy transform.
Use
scipy.stats.levyfor evaluation, simulation, and MLE fitting.See practical workflows: scale inference, hypothesis tests, and a simple Lévy-flight generator.
import platform
import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
import scipy
from scipy import optimize
from scipy.special import erfc, erfcinv
from scipy.stats import chi2, levy, norm
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(7)
print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)
Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0
1) Title & classification#
Name:
levyType: continuous distribution
Support: \(x \in (\mu, \infty)\)
Parameter space: location \(\mu \in \mathbb{R}\) and scale \(c > 0\)
We write:
SciPy uses the same two parameters under the names loc=\mu and scale=c in scipy.stats.levy.
2) Intuition & motivation#
What it models#
The Lévy distribution is a model for strictly-positive, extremely heavy-tailed random variables (up to a location shift): most observations are moderate, but rare samples can be orders of magnitude larger.
A useful tail summary:
PDF: \(f(x) \propto x^{-3/2}\) as \(x \to \infty\)
Survival: \(\mathbb{P}(X>x) \propto x^{-1/2}\) as \(x \to \infty\)
This tail is so heavy that even the first moment diverges.
Typical real-world use cases#
First-passage times: If \(B_t\) is standard Brownian motion and $\(T_a = \inf\{t>0 : B_t = a\},\)\( then \)T_a\( has a Lévy distribution (with \)c=a^2\( and \)\mu=0$). This makes Lévy a natural model for hitting times and barrier-crossing times.
Waiting times with occasional huge delays: queueing-like settings, extreme latency, and other “long right tail” duration phenomena (when a finite-mean model is not appropriate).
Lévy flights: step lengths with a heavy tail can yield trajectories with rare large jumps (anomalous diffusion).
Relations to other distributions#
Inverse-gamma: If \(X \sim \mathrm{Levy}(\mu,c)\) then \(X-\mu\) is an inverse-gamma distribution with shape \(\alpha=\tfrac{1}{2}\) and scale \(\beta=\tfrac{c}{2}\).
Gamma via reciprocal: \(Y = 1/(X-\mu)\) follows \(\mathrm{Gamma}(\tfrac{1}{2},\ \text{rate}=c/2)\).
Normal transform: If \(Z\sim\mathcal{N}(0,1)\) then \(c/Z^2 \sim \mathrm{Levy}(0,c)\).
Stable laws: Lévy is the one-sided stable distribution with index \(\alpha=\tfrac{1}{2}\) (a stable subordinator).
3) Formal definition#
Let \(X \sim \mathrm{Levy}(\mu,c)\) with \(c>0\).
PDF#
For \(x>\mu\):
and \(f(x;\mu,c)=0\) for \(x\le \mu\).
CDF#
For \(x>\mu\):
and \(F(x;\mu,c)=0\) for \(x\le \mu\).
Using the normal CDF \(\Phi\), we also have the equivalent form
Quantile function (inverse CDF)#
For \(p\in(0,1)\):
We’ll implement PDF/CDF (and the quantile) and cross-check against SciPy.
def levy_pdf(x: np.ndarray, loc: float = 0.0, c: float = 1.0) -> np.ndarray:
x = np.asarray(x, dtype=float)
if c <= 0:
raise ValueError("c (scale) must be > 0")
z = x - loc
out = np.zeros_like(z, dtype=float)
m = z > 0
zp = z[m]
out[m] = np.sqrt(c / (2.0 * np.pi)) * np.exp(-c / (2.0 * zp)) / (zp ** 1.5)
return out
def levy_logpdf(x: np.ndarray, loc: float = 0.0, c: float = 1.0) -> np.ndarray:
x = np.asarray(x, dtype=float)
if c <= 0:
raise ValueError("c (scale) must be > 0")
z = x - loc
out = np.full_like(z, -np.inf, dtype=float)
m = z > 0
zp = z[m]
out[m] = 0.5 * (np.log(c) - np.log(2.0 * np.pi)) - 1.5 * np.log(zp) - c / (2.0 * zp)
return out
def levy_cdf(x: np.ndarray, loc: float = 0.0, c: float = 1.0) -> np.ndarray:
x = np.asarray(x, dtype=float)
if c <= 0:
raise ValueError("c (scale) must be > 0")
z = x - loc
out = np.zeros_like(z, dtype=float)
m = z > 0
out[m] = erfc(np.sqrt(c / (2.0 * z[m])))
return out
def levy_ppf(p: np.ndarray, loc: float = 0.0, c: float = 1.0, eps: float = 1e-12) -> np.ndarray:
p = np.asarray(p, dtype=float)
if c <= 0:
raise ValueError("c (scale) must be > 0")
p = np.clip(p, eps, 1.0 - eps)
return loc + c / (2.0 * (erfcinv(p) ** 2))
def levy_ppf_via_norm(p: np.ndarray, loc: float = 0.0, c: float = 1.0, eps: float = 1e-12) -> np.ndarray:
p = np.asarray(p, dtype=float)
if c <= 0:
raise ValueError("c (scale) must be > 0")
p = np.clip(p, eps, 1.0 - eps)
y = norm.ppf(1.0 - p / 2.0) # y>0
return loc + c / (y * y)
# Sanity checks against SciPy
loc, c = 0.2, 1.7
x = loc + np.logspace(-3, 2, 9) * c
p = np.linspace(0.05, 0.95, 7)
print("max |pdf - scipy|:", np.max(np.abs(levy_pdf(x, loc=loc, c=c) - levy.pdf(x, loc=loc, scale=c))))
print("max |cdf - scipy|:", np.max(np.abs(levy_cdf(x, loc=loc, c=c) - levy.cdf(x, loc=loc, scale=c))))
print("max |ppf(erfcinv) - scipy|:", np.max(np.abs(levy_ppf(p, loc=loc, c=c) - levy.ppf(p, loc=loc, scale=c))))
print("max |ppf(norm) - scipy|:", np.max(np.abs(levy_ppf_via_norm(p, loc=loc, c=c) - levy.ppf(p, loc=loc, scale=c))))
max |pdf - scipy|: 5.551115123125783e-17
max |cdf - scipy|: 1.1102230246251565e-16
max |ppf(erfcinv) - scipy|: 1.1368683772161603e-13
max |ppf(norm) - scipy|: 3.552713678800501e-15
4) Moments & properties#
Mean, variance, skewness, kurtosis#
For the Lévy distribution, all ordinary moments of order \(\ge 1/2\) diverge.
In particular:
Mean: \(\mathbb{E}[X]=\infty\) (does not exist as a finite number)
Variance: \(\mathrm{Var}(X)=\infty\)
Skewness / kurtosis: undefined (they require finite variance)
Robust summaries that do exist and are often more meaningful:
Median / quantiles (via the quantile function)
Mode: \(\mu + c/3\)
MGF / Laplace transform and characteristic function#
Because the right tail is so heavy, the MGF does not exist for any \(t>0\).
For \(t<0\), the MGF is a Laplace transform and has a simple closed form:
The characteristic function exists for all real \(t\):
Entropy#
\(X-\mu\) is inverse-gamma with shape \(\alpha=1/2\) and scale \(\beta=c/2\). The (differential) entropy is finite and can be written using the digamma function \(\psi\):
It does not depend on \(\mu\) (shifts do not change differential entropy).
# Visualize the characteristic function φ(t) = exp(i μ t - sqrt(-2 i c t))
mu, c = 0.0, 1.0
t = np.linspace(-25, 25, 3000)
phi = np.exp(1j * mu * t - np.sqrt(-2j * c * t))
fig = make_subplots(rows=1, cols=2, subplot_titles=("Re φ(t)", "Im φ(t)"))
fig.add_trace(go.Scatter(x=t, y=np.real(phi), mode="lines"), row=1, col=1)
fig.add_trace(go.Scatter(x=t, y=np.imag(phi), mode="lines"), row=1, col=2)
fig.update_xaxes(title_text="t", row=1, col=1)
fig.update_xaxes(title_text="t", row=1, col=2)
fig.update_layout(width=950, height=350, showlegend=False)
fig.show()
5) Parameter interpretation#
\(\mu\) (location): shifts the support. Increasing \(\mu\) moves the entire distribution right; \(X-\mu\) is always nonnegative.
\(c\) (scale): controls the typical magnitude and tail heaviness.
Useful facts:
Mode is at \(\mu + c/3\).
Scaling: if \(X \sim \mathrm{Levy}(\mu,c)\) and \(a>0\), then $\(a(X-\mu) \sim \mathrm{Levy}(0,ac),\)\( so quantiles scale linearly with \)c$.
We’ll visualize how changing \(c\) changes the PDF and CDF.
loc = 0.0
c_values = [0.3, 1.0, 3.0]
# Plot as a function of δ = x - loc on a log-x axis to show the heavy tail.
delta = np.logspace(-3, 2, 1200) # 1e-3 ... 1e2
fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))
for c in c_values:
x = loc + c * delta
fig.add_trace(go.Scatter(x=x - loc, y=levy_pdf(x, loc=loc, c=c), mode="lines", name=f"c={c}"), row=1, col=1)
fig.add_trace(go.Scatter(x=x - loc, y=levy_cdf(x, loc=loc, c=c), mode="lines", showlegend=False), row=1, col=2)
fig.update_xaxes(title_text="δ = x - μ", type="log", row=1, col=1)
fig.update_xaxes(title_text="δ = x - μ", type="log", row=1, col=2)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(width=980, height=380)
fig.show()
# The tail is polynomial: f(x) ~ const * (x-μ)^(-3/2)
c = 1.0
loc = 0.0
delta = np.logspace(0, 6, 1200) # 1 ... 1e6
x = loc + delta
pdf = levy_pdf(x, loc=loc, c=c)
# Tail approximation for large δ
pdf_tail = np.sqrt(c / (2.0 * np.pi)) * delta ** (-1.5)
fig = go.Figure()
fig.add_trace(go.Scatter(x=delta, y=pdf, mode="lines", name="pdf"))
fig.add_trace(go.Scatter(x=delta, y=pdf_tail, mode="lines", name=r"√(c/(2π)) δ^{-3/2}", line=dict(dash="dash")))
fig.update_xaxes(title_text="δ = x - μ", type="log")
fig.update_yaxes(title_text="f(x)", type="log")
fig.update_layout(title="Log-log tail: slope -3/2", width=850, height=420)
fig.show()
6) Derivations#
Expectation (why the mean is infinite)#
Work with \(Y=X-\mu \ge 0\). For \(y>0\):
For \(y\ge c\) we have \(\exp\left(-\frac{c}{2y}\right) \ge e^{-1/2}\), so
Then
So \(\mathbb{E}[X]=\mu+\mathbb{E}[Y]=\infty\).
Variance (why it is infinite)#
Similarly,
so the variance is infinite.
Likelihood and MLE for the scale (known location)#
For i.i.d. data \(x_1,\dots,x_n\) with \(x_i>\mu\), the log-likelihood is
If \(\mu\) is known, differentiate w.r.t. \(c\) and set to zero:
Estimating \(\mu\) by MLE is tricky: with enough samples the likelihood typically increases as \(\mu \uparrow \min_i x_i\) (a common issue with support parameters).
def levy_loglik(x: np.ndarray, loc: float, c: float) -> float:
x = np.asarray(x, dtype=float)
return float(np.sum(levy_logpdf(x, loc=loc, c=c)))
def levy_mle_scale_given_loc(x: np.ndarray, loc: float) -> float:
x = np.asarray(x, dtype=float)
z = x - loc
if np.any(z <= 0):
raise ValueError("All observations must satisfy x_i > loc")
return float(x.size / np.sum(1.0 / z))
# Demonstration: recover c when loc is known
true_loc, true_c = 0.0, 1.5
x = levy.rvs(loc=true_loc, scale=true_c, size=2000, random_state=rng)
c_hat = levy_mle_scale_given_loc(x, loc=true_loc)
print("true c:", true_c)
print("MLE c (loc known):", c_hat)
# Compare to SciPy's fit with fixed loc
loc_fit, c_fit = levy.fit(x, floc=true_loc)
print("SciPy fit (loc fixed):", (loc_fit, c_fit))
true c: 1.5
MLE c (loc known): 1.5289999728329602
SciPy fit (loc fixed): (0.0, 1.5290039062500016)
# Profile likelihood over loc (illustration only)
def levy_profile_nll_loc(x: np.ndarray, loc: float) -> float:
x = np.asarray(x, dtype=float)
if loc >= np.min(x):
return np.inf
try:
c_hat = levy_mle_scale_given_loc(x, loc=loc)
except ValueError:
return np.inf
return -levy_loglik(x, loc=loc, c=c_hat)
x = levy.rvs(loc=0.3, scale=1.0, size=200, random_state=rng)
min_x = float(np.min(x))
# Search in a window below the minimum observation
upper = min_x - 1e-6
lower = min_x - 5.0
res = optimize.minimize_scalar(
lambda loc: levy_profile_nll_loc(x, loc=loc),
bounds=(lower, upper),
method="bounded",
)
loc_hat = float(res.x)
c_hat = levy_mle_scale_given_loc(x, loc=loc_hat)
print("min(x):", min_x)
print("profile-MLE loc:", loc_hat)
print("profile-MLE c:", c_hat)
print("SciPy fit:", levy.fit(x))
min(x): 0.41549549648756096
profile-MLE loc: 0.30890712001029735
profile-MLE c: 1.0516969535184806
SciPy fit: (0.3089117456740079, 1.0516693766560499)
7) Sampling & simulation (NumPy-only)#
A convenient exact sampler comes from the identity:
So to sample \(X \sim \mathrm{Levy}(\mu,c)\):
Draw \(Z \sim \mathcal{N}(0,1)\)
Return \(X = \mu + c/Z^2\)
This uses only NumPy and is typically faster than inverse-CDF sampling (and avoids erfcinv).
def sample_levy_numpy(
rng: np.random.Generator,
size: int,
loc: float = 0.0,
c: float = 1.0,
) -> np.ndarray:
if c <= 0:
raise ValueError("c (scale) must be > 0")
z = rng.standard_normal(size)
# Extremely unlikely but possible with finite-precision RNG outputs.
# Resample any exact zeros to avoid division by zero.
while True:
m = z == 0.0
if not np.any(m):
break
z[m] = rng.standard_normal(int(np.sum(m)))
return loc + c / (z * z)
# Compare NumPy-only sampler to SciPy on quantiles
n = 200_000
loc, c = 0.0, 1.0
x_np = sample_levy_numpy(rng, n, loc=loc, c=c)
qs = [0.1, 0.25, 0.5, 0.75, 0.9]
q_emp = np.quantile(x_np, qs)
q_the = levy.ppf(qs, loc=loc, scale=c)
print("quantiles:")
for p, qe, qt in zip(qs, q_emp, q_the):
print(f" p={p:>4}: empirical={qe:>10.4f} theory={qt:>10.4f}")
quantiles:
p= 0.1: empirical= 0.3699 theory= 0.3696
p=0.25: empirical= 0.7591 theory= 0.7557
p= 0.5: empirical= 2.2115 theory= 2.1981
p=0.75: empirical= 9.9013 theory= 9.8492
p= 0.9: empirical= 63.8245 theory= 63.3281
8) Visualization#
We’ll visualize:
the PDF and CDF
a Monte Carlo histogram with PDF overlay
the running mean (which does not stabilize because the mean is infinite)
loc, c = 0.0, 1.0
# Use log-spaced δ to cover both the peak and far tail.
delta = np.logspace(-3, 3, 2000)
x = loc + delta
pdf = levy_pdf(x, loc=loc, c=c)
cdf = levy_cdf(x, loc=loc, c=c)
fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))
fig.add_trace(go.Scatter(x=delta, y=pdf, mode="lines", name="pdf"), row=1, col=1)
fig.add_trace(go.Scatter(x=delta, y=cdf, mode="lines", name="cdf", showlegend=False), row=1, col=2)
fig.update_xaxes(title_text="δ = x - μ", type="log", row=1, col=1)
fig.update_xaxes(title_text="δ = x - μ", type="log", row=1, col=2)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(width=980, height=380)
fig.show()
# Monte Carlo histogram (use a high quantile clip to keep the plot readable)
n = 200_000
x = sample_levy_numpy(rng, n, loc=0.0, c=1.0)
clip_q = 0.995
x_clip = np.quantile(x, clip_q)
x_plot = x[x <= x_clip]
# Histogram on log-x scale
bins = np.logspace(-3, np.log10(x_clip), 80)
hist, edges = np.histogram(x_plot, bins=bins, density=True)
centers = np.sqrt(edges[:-1] * edges[1:])
pdf_centers = levy_pdf(centers, loc=0.0, c=1.0)
fig = go.Figure()
fig.add_trace(go.Bar(x=centers, y=hist, name="MC histogram", marker=dict(opacity=0.55)))
fig.add_trace(go.Scatter(x=centers, y=pdf_centers, mode="lines", name="theoretical pdf"))
fig.update_xaxes(title_text="x (log scale)", type="log")
fig.update_yaxes(title_text="density")
fig.update_layout(
title=f"Histogram (clipped at {clip_q:.3f} quantile to {x_clip:.2f})",
width=900,
height=450,
)
fig.show()
# Running mean does not stabilize (mean is infinite)
n = 60_000
x = sample_levy_numpy(rng, n, loc=0.0, c=1.0)
run_mean = np.cumsum(x) / np.arange(1, n + 1)
fig = go.Figure()
fig.add_trace(go.Scatter(x=np.arange(1, n + 1), y=run_mean, mode="lines"))
fig.update_xaxes(title_text="n", type="log")
fig.update_yaxes(title_text="running mean", type="log")
fig.update_layout(title="Running mean keeps drifting upward", width=900, height=420)
fig.show()
9) SciPy integration (scipy.stats.levy)#
SciPy provides the Lévy distribution as scipy.stats.levy.
Common methods:
levy.pdf(x, loc, scale)levy.cdf(x, loc, scale)levy.rvs(loc, scale, size, random_state)levy.fit(data)(MLE forlocandscale)
We’ll do a small end-to-end example.
# End-to-end: simulate, fit, and compare
true_loc, true_c = 0.2, 1.0
x = levy.rvs(loc=true_loc, scale=true_c, size=5000, random_state=rng)
loc_hat, c_hat = levy.fit(x)
loc_hat_fixed, c_hat_fixed = levy.fit(x, floc=true_loc)
print("true (loc, c):", (true_loc, true_c))
print("fit (loc, c):", (loc_hat, c_hat))
print("fit with loc fixed:", (loc_hat_fixed, c_hat_fixed))
# Compare a few values
grid = true_loc + np.logspace(-3, 2, 8)
print("pdf(grid) from our implementation:", levy_pdf(grid, loc=true_loc, c=true_c))
print("pdf(grid) from SciPy:", levy.pdf(grid, loc=true_loc, scale=true_c))
true (loc, c): (0.2, 1.0)
fit (loc, c): (0.20295360934385015, 0.9358360873038567)
fit with loc fixed: (0.2, 0.9446289062499998)
pdf(grid) from our implementation: [0. 0. 0. 0.2108 0.3262 0.0485 0.0046 0.0004]
pdf(grid) from SciPy: [0. 0. 0. 0.2108 0.3262 0.0485 0.0046 0.0004]
10) Statistical use cases#
Hypothesis testing: exact inference for the scale (known location)#
If the location \(\mu\) is known, define
Then \(Y_i \sim \mathrm{Gamma}(\tfrac{1}{2},\ \text{rate}=c/2)\) and the sum
satisfies an exact chi-square identity:
So you can test \(H_0:c=c_0\) or build an exact confidence interval for \(c\).
Bayesian modeling: conjugate update for the scale (known location)#
With known \(\mu\), the likelihood in \(c\) is
so a Gamma prior on \(c\) is conjugate:
Generative modeling: Lévy flights#
Use Lévy-distributed step lengths to generate paths with rare large jumps (a toy model for anomalous diffusion).
# Exact confidence interval / test for c when loc is known
def levy_scale_ci_known_loc(x: np.ndarray, loc: float, alpha: float = 0.05) -> tuple[float, float]:
x = np.asarray(x, dtype=float)
z = x - loc
if np.any(z <= 0):
raise ValueError("All observations must satisfy x_i > loc")
S = float(np.sum(1.0 / z))
n = x.size
lo = chi2.ppf(alpha / 2.0, df=n) / S
hi = chi2.ppf(1.0 - alpha / 2.0, df=n) / S
return float(lo), float(hi)
true_loc, true_c = 0.0, 1.5
x = levy.rvs(loc=true_loc, scale=true_c, size=400, random_state=rng)
ci = levy_scale_ci_known_loc(x, loc=true_loc, alpha=0.05)
c_hat = levy_mle_scale_given_loc(x, loc=true_loc)
print("true c:", true_c)
print("MLE c:", c_hat)
print("95% CI for c (exact):", ci)
print("contains true c?", ci[0] <= true_c <= ci[1])
true c: 1.5
MLE c: 1.4432416413976406
95% CI for c (exact): (1.2501422793918129, 1.650005786032107)
contains true c? True
# Bayesian posterior over c with a Gamma prior (known loc)
from scipy.stats import gamma
true_loc, true_c = 0.0, 1.0
x = levy.rvs(loc=true_loc, scale=true_c, size=200, random_state=rng)
z = x - true_loc
S = float(np.sum(1.0 / z))
n = x.size
# Prior: c ~ Gamma(a0, rate=b0)
a0, b0 = 2.0, 1.0
# Posterior: c | x ~ Gamma(a_post, rate=b_post)
a_post = a0 + n / 2.0
b_post = b0 + 0.5 * S
# SciPy uses 'scale' = 1/rate
prior = gamma(a=a0, scale=1.0 / b0)
post = gamma(a=a_post, scale=1.0 / b_post)
grid = np.linspace(0.001, post.ppf(0.995), 1200)
fig = go.Figure()
fig.add_trace(go.Scatter(x=grid, y=prior.pdf(grid), mode="lines", name="prior"))
fig.add_trace(go.Scatter(x=grid, y=post.pdf(grid), mode="lines", name="posterior"))
fig.add_vline(x=true_c, line=dict(dash="dash"), annotation_text="true c")
fig.update_xaxes(title_text="c")
fig.update_yaxes(title_text="density")
fig.update_layout(title="Posterior over scale c (known loc)", width=900, height=420)
fig.show()
print("posterior mean c:", post.mean())
print("posterior sd c:", post.std())
posterior mean c: 1.0231191657045369
posterior sd c: 0.1013038928094692
# Lévy flight: 2D random walk with Lévy step lengths
n_steps = 1500
c = 0.02 # smaller c -> smaller typical steps (but still heavy-tailed)
lengths = sample_levy_numpy(rng, n_steps, loc=0.0, c=c)
angles = rng.uniform(0.0, 2.0 * np.pi, size=n_steps)
dx = lengths * np.cos(angles)
dy = lengths * np.sin(angles)
x = np.cumsum(dx)
y = np.cumsum(dy)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode="lines", line=dict(width=1)))
fig.update_xaxes(title_text="x")
fig.update_yaxes(title_text="y", scaleanchor="x", scaleratio=1)
fig.update_layout(title="Toy Lévy flight (occasional big jumps)", width=850, height=650)
fig.show()
11) Pitfalls#
Invalid parameters: require \(c>0\). The distribution is undefined for nonpositive scale.
Support violations: densities and likelihoods assume \(x>\mu\); if you fit/optimize \(\mu\), any \(x_i\le\mu\) makes the likelihood \(-\infty\).
Infinite moments: sample mean/variance are unstable; moment-matching is not appropriate.
Visualization needs clipping or log axes: a few extreme draws can hide the bulk of the distribution.
Sampling can create huge values: \(X=\mu + c/Z^2\) explodes when \(|Z|\) is very small (this is real, not a bug). Use float64 and be careful with downstream computations.
Fitting
loccan stick to the sample minimum: MLEs for support parameters often land on (or near) the boundary; consider fixinglocif it is known from the problem.Prefer
logpdffor likelihoods: products of PDFs underflow; sums of log-PDFs are stable.
12) Summary#
levyis a continuous, one-sided heavy-tailed distribution on \((\mu,\infty)\) with parameters \((\mu,c)\).It has closed-form PDF/CDF/quantile (involving
erfc) and a simple NumPy-only sampler: \(X=\mu+c/Z^2\).The tail is polynomial (\(\sim x^{-3/2}\)), implying infinite mean and variance.
It is closely related to inverse-gamma and to Brownian first-passage times.
With known \(\mu\), the scale \(c\) admits exact inference via the identity \(c\sum 1/(X_i-\mu)\sim\chi^2_n\).
References:
SciPy docs:
scipy.stats.levyFeller — An Introduction to Probability Theory and Its Applications (first-passage times)
Sato — Lévy Processes and Infinitely Divisible Distributions